home *** CD-ROM | disk | FTP | other *** search
/ Dr. Windows 3 / dr win3.zip / dr win3 / BUSINESS / XLMATH22.ZIP / SOURCE.ZIP / XLMATH.C < prev    next >
C/C++ Source or Header  |  1993-04-19  |  21KB  |  851 lines
  1. /*----------------------------------------------------------------------*\
  2.  
  3.  XLMATH    
  4.  ========
  5.  Copyright    1992 by Roy Kari
  6.  
  7.  Author:    Roy Kari
  8.              Dept. of Chemistry
  9.             Laurentian University
  10.             Sudbury, Ont.
  11.             Canada        P3E 2C6
  12.             (705) 675-1151
  13.             Internet: "ROY@NICKEL.LAURENTIAN.CA"
  14.  
  15.  Revision history:
  16.  
  17.     1.0    September, 1992
  18.         initial version
  19.  
  20.     2.0    January, 1993
  21.         added dialog boxes
  22.  
  23.     2.1    April, 1993
  24.         bug fixes
  25.  
  26.  
  27.     This module contains routines written to interface normal C
  28.     routines to an Excel custom function format. All custom
  29.     functions are assumed to pass a floating point array to the custom
  30.     function, and the custom function returns an XLOPER, usually an
  31.     array. If you wish to use a dialog box, then the dialog routine
  32.     is written to interface the custom function to a dialog box.
  33.     The dialog interface must get the input from the Excel sheet,
  34.     change it to a floating point array, and then call the custom function.
  35.     On return from the custom function, the dialog interface must
  36.     display the return XLOPER in the sheet, and free all allocated
  37.     memory.
  38.  
  39.     2nd thoughts;
  40.     If I had this to over again, I would make the custom functions
  41.     accept their input in the form of XLOPER's rather than a 
  42.     floating point array. This would eliminate the need to convert.
  43.  
  44. \*----------------------------------------------------------------------*/
  45.  
  46.  
  47. /* --------------------------< Include files >--------------------- */
  48.  
  49. #include <windows.h>
  50. #include <xlcall.h>
  51.  
  52. #include <float.h>
  53. #include <math.h>
  54. #include <stdlib.h>
  55.  
  56. #include <framewrk.h>
  57. #include "xlmath.h"
  58. #include "xlutil.h"
  59. #include "xlcurve.h"
  60. #include "xlmcfit.h"
  61.  
  62. REALTYPE predict (int k, PARVEC p, REALTYPE xk);
  63.  
  64. /********************************************************************
  65.  Diagonalize()
  66.  =============
  67.  Diagonalize a real symmetric matrix. Eigenvalues are returned
  68.  *******************************************************************/
  69. LPXLOPER PASCAL FAR __export Diagonalize(LPFP lpHmat)
  70. {
  71.     LPMULTI lpMulti = NULL;    // return XLOPER's
  72.     LPLPREAL lplpHmat;        // pointers for EXCEL allocated array
  73.     LPREAL lpEivec;            // eigenvector array
  74.     LPLPREAL lplpEivec;        // pointers for above
  75.     LPREAL lpEival;            // eigenvalue array
  76.  
  77.     WORD i,j,wRowIndex;        // index
  78.     int nIter;                // no iterations in _Diagonalize()
  79.     WORD wRows, wCols;        // saves typing
  80.     BOOL Success = TRUE;
  81.  
  82.     wRows = lpHmat->wRows;
  83.     wCols = lpHmat->wCols;
  84.  
  85.     // error if not square array
  86.     if (wRows != wCols)
  87.     {
  88.         ErrorHandler(XLM_NOT_SQUARE);
  89.         return (glpxlError);
  90.     }
  91.     
  92.     //     allocate a Global 2-D array for eigenvectors. This array allocation
  93.     //    demonstrates the use of 2-D array indexing M[i][j]. It is 
  94.     //    required in this case because the working routine _Diagonalize() was
  95.     //    written in this manner and it was too much work to convert
  96.     //    it to anything else.
  97.     if ((lpEivec = (LPREAL)GetMem( sizeof(double)*wRows*wCols)) == NULL)
  98.     {
  99.         return (glpxlError);
  100.  
  101.     }
  102.     if ((lplpEivec = InitPointers(lpEivec, wRows, wCols)) == NULL)
  103.     {
  104.         Success = FALSE;
  105.         goto Error3;        // free lpEivec
  106.     }
  107.  
  108.     // create pointers for Hmat; array allocated in Excel
  109.     if ((lplpHmat = InitPointers(&lpHmat->Data[0], wRows, wRows))
  110.             == NULL)
  111.  
  112.     {
  113.         Success = FALSE;
  114.         goto Error2;        // free lpEivec & lplpEivec
  115.     }
  116.  
  117.     // allocate a Global 1-D array for eigenvalues
  118.     if ((lpEival = (LPREAL)GetMem( sizeof(double)*wRows )) == NULL)
  119.     {
  120.         Success = FALSE;
  121.         goto Error1;        // free lpEivec, lplpEivec & lplpHmat
  122.     }
  123.  
  124.     // diagonalize the matrix
  125.     nIter = _Diagonalize(lplpHmat, lplpEivec, lpEival, wRows);
  126.  
  127.  
  128.     //
  129.     // return vales to EXCEL
  130.     //
  131.     
  132.     // init MULTI 1 row larger to store eigenvalues in last row
  133.     if ((lpMulti =InitMulti(wRows+1, wRows, xltypeNum)) == NULL)
  134.     {
  135.         Success = FALSE;
  136.     }
  137.     else
  138.     {
  139.         // copy the eigenvectors
  140.         for (i = 0; i < wRows; i++)
  141.         {
  142.             wRowIndex = i*wCols;
  143.             for (j=0; j < wCols; j++)
  144.             {
  145.                 lpMulti->Data[wRowIndex+j].val.num = lplpEivec[i][j];
  146.             }
  147.         }
  148.  
  149.         // copy the eigenvalues
  150.         wRowIndex = wRows*wCols;
  151.         for (j = 0; j < wCols; j++)
  152.         {
  153.             lpMulti->Data[wRowIndex+j].val.num = lpEival[j];
  154.         }
  155.     }
  156.  
  157.     // free pointers and arrays
  158.     FreeMem((LPVOID)lpEival);
  159. Error1:
  160.     MemFreePtr((LPVOID)lplpHmat);
  161. Error2:
  162.     MemFreePtr((LPVOID)lplpEivec);
  163. Error3:
  164.     FreeMem((LPVOID)lpEivec);
  165.     if (!Success)
  166.     {
  167.         return (glpxlError);
  168.     }
  169.  
  170.     // return XLOPER to EXCEL
  171.     return ( (LPXLOPER)lpMulti);
  172. }
  173.  
  174. /********************************************************************
  175.  MODensity()
  176.  =============
  177.  Compute the MO density as C(T)xOxC
  178.  *******************************************************************/
  179. LPXLOPER FAR PASCAL __export MODensity(LPFP lpEivec, LPFP lpOcc)
  180. {
  181.     LPMULTI lpDen = NULL;
  182.     UINT i, j, k;
  183.     UINT i_Idx, j_Idx, ij_Idx, NumVec;
  184.     BOOL bError = FALSE;
  185.  
  186.     NumVec = lpEivec->wRows;
  187.  
  188.     if (lpOcc->wRows == 1)
  189.     {
  190.         // Occ is row vector (1xN)
  191.         if (lpOcc->wCols != NumVec)
  192.         {
  193.             bError = TRUE;
  194.         }
  195.     }
  196.     else
  197.     {
  198.         // Occ is a column vector (Nx1)
  199.         if ( (lpOcc->wRows != NumVec) || (lpOcc->wCols != 1) )
  200.         {
  201.             bError = TRUE;
  202.         }
  203.     }
  204.     if (bError)
  205.     {
  206.         ErrorHandler(XLM_MODENSITY);
  207.         return (glpxlError);
  208.     }
  209.  
  210.     if ((lpDen =InitMulti(NumVec, NumVec, xltypeNum)) == NULL)
  211.     {
  212.         return (glpxlError);
  213.     }
  214.  
  215.     i_Idx = 0;
  216.  
  217.     for (i = 0; i < NumVec; i++)
  218.     {
  219.         j_Idx = 0;
  220.         for (j = 0; j < NumVec; j++)
  221.         {
  222.             ij_Idx = i_Idx+j;
  223.             lpDen->Data[ij_Idx].val.num = 0.0;
  224.             for (k = 0; k < NumVec; k++)
  225.             {
  226.                 // den[i][j] = Eivec[i][k]*Eivec[j][k]*Occ[k]
  227.                 lpDen->Data[ij_Idx].val.num += lpEivec->Data[i_Idx+k] *
  228.                                                 lpEivec->Data[j_Idx+k] *
  229.                                                 lpOcc->Data[k];
  230.             }
  231.             j_Idx += NumVec;
  232.         }
  233.         i_Idx += NumVec;
  234.     }
  235.     return ((LPXLOPER)lpDen);
  236. }
  237.  
  238.  
  239. /********************************************************************
  240.  PolyCurveFit()
  241.  =============
  242.  This function is a generalized least squares curve fitting function.
  243.  It fits a polynomial with linear coefficients to a dependent -
  244.  independent variable set of data
  245.  *******************************************************************/
  246. LPXLOPER PASCAL FAR _export PolyCurveFit(LPFP lpIndVar, LPFP lpDepVar,
  247.                 unsigned int Order)
  248. {
  249. #define    nCols    3    // placed in a NumObs x nCols array
  250.  
  251.     LPMULTI lpMulti;
  252.     LPREAL     lpYest, lpResid, lpPolycoef, lpCoefsig;
  253.     static double    See, Rsqrval, Rval, dferror;
  254.     static WORD wNumObs;
  255.     WORD wRows = lpIndVar->wRows;
  256.     WORD wCols = lpIndVar->wCols;
  257.     WORD wMinRows = 2*(Order+1) + 4;    // needed for return values
  258.     BOOL Success = TRUE;
  259.     BOOL bColumnVector = TRUE;
  260.     WORD i, wRow1, wRow2, wRowIndex;
  261.  
  262.     wNumObs = wRows;
  263.     if (wCols > wRows)
  264.     {
  265.         // change to row vectors
  266.         bColumnVector = FALSE;
  267.         wNumObs = wCols;
  268.     }
  269.  
  270.     // check input dimensions so both same
  271.     if     (     (lpDepVar->wRows != lpIndVar->wRows)     ||
  272.             (lpDepVar->wCols != lpIndVar->wCols)    ||
  273.             (Order >= wNumObs - 1)
  274.         )
  275.     {
  276.         ErrorHandler(XLM_POLY);
  277.         return(glpxlError);
  278.     }
  279.  
  280.     if (wMinRows < wNumObs)
  281.         wMinRows = wNumObs;
  282.  
  283.  
  284.     // memory for estimated Y values
  285.     lpYest = (LPREAL)GetMem( wNumObs*sizeof(double));
  286.     if (lpYest == NULL)
  287.     {
  288.         Success = FALSE;
  289.         goto Error4;
  290.     }
  291.     // memory for residuals
  292.     if ((lpResid = (LPREAL)GetMem( wNumObs*sizeof(double))) == NULL)
  293.     {
  294.         Success = FALSE;
  295.         goto Error3;
  296.     }
  297.  
  298.     //memory for coefficients of fitted polynomial
  299.     if ((lpPolycoef = (LPREAL)GetMem( (Order+1)*sizeof(double))) == NULL)
  300.     {
  301.         Success = FALSE;
  302.         goto Error2;
  303.     }
  304.  
  305.     // memory for standard errors of coefficient estimates
  306.     if ((lpCoefsig = (LPREAL)GetMem( (Order+1)*sizeof(double))) == NULL)
  307.     {
  308.         Success = FALSE;
  309.         goto Error1;
  310.     }
  311.  
  312.     Success  = _PolyCurveFit((LPREAL)&lpIndVar->Data[0],
  313.         (LPREAL)&lpDepVar->Data[0], (int)wNumObs, (int)Order, lpPolycoef,
  314.         lpYest, lpResid, (NPREAL)&See, lpCoefsig, &Rsqrval, &Rval, &dferror);
  315.     if (!Success)
  316.         goto Error0;
  317.  
  318.     // init XLOPER
  319.     if (bColumnVector)
  320.         lpMulti = InitMulti(wMinRows, nCols, xltypeNum);
  321.     else
  322.         lpMulti = InitMulti(nCols, wMinRows, xltypeNum);
  323.     if (lpMulti==NULL)
  324.     {
  325.         Success = FALSE;
  326.         goto Error0;
  327.     }
  328.  
  329.     // stuff return data into Multi
  330.     if (bColumnVector)
  331.